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Abstract. A set-up is introduced which can be super- 
imposed onto the existing solar flare cellular automata 
(CA) models, and which specifies the interpretation of the 
model's variables. It extends the CA models, yielding the 
magnetic field, the current, and an approximation to the 
electric field, in a way that is consistent with Maxwell's 
and the MHD equations. Applications to several solar flare 
CA models during their natural state (self-organized crit- 
icality (SOC)) show, among others, that (1) the magnetic 
field exhibits characteristic large-scale organization over 
the entire modeled volume; (2) the magnitude of the cur- 
rent seems spatially dis-organized, with no obvious ten- 
dency towards large-scale structures or even local organi- 
zation; (3) bursts occur at sites with increased current, 
and after a burst the current is relaxed; (4) by estimating 
the energy released in individual bursts with the use of the 
current as Ohmic dissipation, it turns out that the power- 
law distributions of the released energy persist. The CA 
models, extended with the set-up, can thus be considered 
as models for energy-release through current- dissipation. 
The concepts of power-law loading and anisotropic events 
(bursts) in CA models are generalized to 3-D vector-field 
models, and their effect on the magnetic field topology is 
demonstrated. 

Key words: solar flares, cellular automata, MHD, non- 
linear processes, chaos, turbulence 



1. Introduction 

Cellular automata (CA) models for solar flares are suc- 
cessful in explaining solar flare statistics (peak flux, to- 
tal flux, and duration distributions; Lu & Hamilton 1991 
(hereafter LH91); Lu et al. 1993; Vlahos et al. 1995; Geor- 



goulis & Vlahos 1996, 1998; Galsgaard 1996). They sim- 
plify strongly the details of the involved physical pro- 
cesses, and achieve in this way to model large volumes with 
complex field topologies and a large number of events. On 
the other hand, MHD simulations give insight into the de- 
tails of the local processes, they are limited, however, to 
modeling relatively small fractions of active regions, due to 
the lack of computing power, yielding thus poor statistics 
and difficulties in comparing results to observations (e.g. 
Mikic et al. 1989; Strauss 1993; Longcope & Sudan 1994; 
Einaudi et al. 1996; Galsgaard & Nordlund 1996; Hendrix 
& Van Hoven 1996; Dmitruk & Gomez 1998; Galtier & 
Pouquet 1998; Georgoulis et al. 1998; Karpen et al. 1998; 
Einaudi & Velh 1999). The global MHD flare models are 
still in the state of rather qualitative flare scenarios. 

The MHD and the CA approach to solar flares seem to 
have very little in common: The former are a set of partial 
differential equations, based on fluid-theory and Maxwell's 
equations, whereas the latter are a set of abstract evo- 
lution rules, based (in the case of solar flares) on the 
analogy to critical phenomena in (theoretical) sand-piles. 
The scope of this paper is to bridge the gap in-between 
these two approaches: the solar flare CA models are re- 
interpreted and extended so as (i) to make these mod- 
els completely compatible with MHD and with Maxwell's 
equations, and so that (ii) all relevant MHD variables 
are made available (e.g. the current and the electric field, 
which so far were not available in CA models). 

In an earlier paper (Isliker et al. 1998), we have ana- 
lyzed the existing solar flare CA models for their sound- 
ness with MHD. We asked the question whether the fields 
in these CA models and the evolution rules can be inter- 
preted in terms of MHD. It turned out that these models 
can indeed be interpreted as a particular way of imple- 
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mcnting numerically the MHD equations. This fact is not 
trivial, since these models had been derived in quite close 
analogy to the sand-pile CA model of Bak et al. (1987 and 
1988), with vague association of the model's variables with 
physical quantities. For instance, some authors (Lu et al. 
1993) explicitly discuss the question whether their basic 
grid variable is the magnetic field or not, without reaching 
to a definite conclusion. Isliker et al. (1998) brought forth 
not only how the existing CA models are related to MHD 
and what simplifications arc hidden, but also where they 
differ from or even violate the laws of MHD and Maxwell's 
equation. Important is the fact that though the existing 
CA models can be considered as a strongly simplified nu- 
merical solution of the (simplified) MHD equations, they 
do not represent the discretized MHD equations: the time- 
step and the spacing between two grid sites are not small 
(in a physical sense), but finite; they are a typical tem- 
poral and spatial scale of the diffusive processes involved 
(see Isliker ct al. 1998). 

From the point of view of MHD, the main short- 
comings of the existing CA models are (Isliker et al. 1998): 

(1) There is no control over consistency with Maxwell's 
equations. Interpreting, for instance, the vector-field in 
the CA models as the magnetic field leads to the problem 
that the gradient of the field (VB) cannot be controlled. 

(2) Secondary quantities, such as currents, are not avail- 
able, and they cannot be introduced in the straightforward 
way by replacing differential expressions by difference- 
expressions, since, as mentioned, the grid-size must be 
considered finite (see also App. B.l). This lack of know- 
ing how to calculate derivatives made it also useless to 
interpret the primary vector-field in the CA models as 
the vector potential (to avoid the VB-problem), since B 
could not be derived. The physical interpretation of these 
CA models remained so far problematic. 

There are two basically different ways of developing 
CA models for flares further: (i) Either one considers CA 
models per se, tries to change the existing models further 
or invent new ones, with the only aim of adjusting them 
to reproduce still better the observations, i.e. one makes 
them a tool the results of which explain and maybe pre- 
dict observed properties of flares. In this approach, one 
has not to care about possible inconsistencies with MHD 
or even Maxwell's equations, the various components of 
the model are purely instrumentalistic. (ii) On the other 
hand, one may care about the physical identification and 
interpretation of the various components of the model, not 
just of its results, and one may want the CA model to be- 
come consistent with the other approach to solar flares, 
namely MHD. In the approach (ii), some of the freedom 
one has in constructing CA models will possibly be re- 
duced, since there are more 'boundary conditions' to be 
fulfilled in the construction of the model: the observations 
must be reproduced and consistency with MHD has to be 
reached. (Trials to construct new CA models which are 
based on MHD and not on the sand-pile analogy were re- 



cently made by Einaudi & Velli (1999), MacPhcrson & 
MacKinnon (1999), Longcope and Noonan (2000), and Is- 
liker et al. (2000a).) 

Our aim is in-between these two alternatives: we con- 
struct a set-up which can be superimposed onto each clas- 
sical solar flare CA model, and which makes the latter 
intcrpretable in a MHD-consistcnt way (by classical CA 
models we mean the models of LH91, Lu et al. 1993, Vla- 
hos et al. 1995, Georgoulis & Vlahos 1996, 1998, Galsgaard 
1996, and their modifications, which are based on the 
sand-pile analogy). The set-up thus specifies the physical 
interpretation of the grid-variables and allows the deriva- 
tion of quantities such as currents etc. It does not interfere 
with the dynamics of the CA (unless wished): loading, re- 
distributing (bursting), and the appearance of avalanches 
and self-organized criticality (SOC), if the latter are im- 
plied by the evolution rules, remain unchanged. The result 
is therefore still a CA model, with all the advantages of 
CA, namely that they are fast, that they model large spa- 
tial regions (and large events), and therewith that they 
yield good statistics. Since the set-up introduces all the 
relevant physical variables into the context of the CA mod- 
els, it automatically leads to a better physical understand- 
ing of the CA models. It reveals which relevant plasma 
processes and in what form are actually implemented, and 
what the global flare scenario is the CA models imply. All 
this was more or less hidden so far in the abstract evo- 
lution rules. It leads also to the possibility to change the 
CA models (the rules) at the guide-line of MHD, if this 
should become desirable. Not least, the set-up opens a way 
for further comparison of the CA models to observations. 

In Sec. 2, we introduce our set-up. Applying it to sev- 
eral CA models (Sec. 3), we will demonstrate the useful- 
ness and some of the beneflts such extended models (i.e. 
classical models extended with our set-up) provide over 
the classical CA models, and we will reveal basic physi- 
cal features of the CA models. The potential of the ex- 
tended models to explain more observational facts than 
the classical CA models is, among others, outlined in the 
conclusions (Sec. 4). 

2. Introduction of the set-up 

The set-up we propose can be superimposed onto solar 
flare CA models which use a 3-D grid and a basic 3~D 
vector grid- variable, say A. The corresponding set of evo- 
lution rules is not changed. (With a few modifications, the 
set-up can also be superimposed onto CA models which 
use a scalar field in a planar grid, which our set-up nec- 
essarily interprets as slab geometry, as will become clear 
later.) 

We introduce our model on the example of the solar 
flare CA model of LH91, which we summarize here in order 
to make the subsequent presentation more concrete: 
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2.1. Summary of the CA model of LH91 

In the LH91 model, to each grid-site Xijk of a 3 D cubic 
grid a 3-D vector A^fe is assigned. Initially, Aij^ is set to 
(1,1, 1)-^, everywhere. The system is then loaded with the 
repeated dropping of increments at randomly chosen sites 

Xijk (one per time-step) 

A{t + 1, ajy-fe) = -4(t, Xijk) + SA{t, Xijk), (1) 

where SA{t, Xijk) has all its components as random num- 
bers uniformly distributed in [—0.2, 0.8]. 

After each loading event, the system is checked for 
whether the local 'stress', defined as 

dAijk '.= Aijk — ^ ^ An.n.j (2) 

n.n. 

where the sum goes over the six nearest neighbours of the 
central point Xijk, exceeds a threshold Acr, i.e. whether 

\dAijk\>Acr, (3) 

where Acr = 7 is used. If this is the case, the field in the 
neighbourhood of the critical site is redistributed accord- 
ing to 

■Aijk ^ A-ijk —dAijk (4) 
for the central point, and 

A-n.n. ^ A.n,n. ~\~ "^dAijk (5) 

for its six nearest neighbours. In such a redistribution 
event (burst), energy amounting to 

E^r = l\dAijk\' (6) 

is assumed to be released. The grid is scanned again and 

again to search for second, third etc. generation bursts, 
until the system is nowhere critical anymore and returns 
to the loading phase (the details we apply concerning the 
temporal evolution of the model arc given in App. A; they 
are not explicitly stated in LH91). The field outside the 
grid is held constant and assumed to be zero. 

2.2. Our set-up 

We turn now to introducing our set-up, starting with a 
specification: We interpret the vector Aijk at the grid sites 
Xijk to denote the local vector-field, A{xijk). Note that 
this was not specified in the classical CA models. Lu ct al. 
(1993) for instance discuss this point: it might also have 
been thought of as a mean local field, i.e. the average over 
an elementary cell in the grid. 

Guided by the idea that we want to assure VJB = 
for the magnetic field B, which is most easily achieved 
by having the vector-potential A as the primary variable 
and letting B be the corresponding derivative of A {B = 
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V A A) , we furthermore assume that the grid variable A 
of the CA model is identical with the vector-potential. 

The remaining and actually most basic problem then 
is to find an adequate way to calculate derivatives in the 
grid. In general, CA models assume that the grid-spacing 
is finite, which also holds for the CA model of LH91 (as 
shown in detail by Isliker et al. 1998), so that the most 
straightforward way of replacing differential expressions 
with difference expressions is not adequate (see the de- 
tailed discussion in App. B.l, below; Vassiliadis ct al. 
(1998) suggested to interpret CA models as the straight- 
forwardly discretized (simplified) MHD equations, which 
we find problematic for the reasons given in App. B.l, and 
we therefore do not follow this approach, here). 

Consequently, one has to find a way of continuing the 
vector-field into the space in-between the grid-sites, which 
will allow to calculate derivatives. There is, of course, an 
infinite number of possibilities to do so, and the problem 
cannot have a unique solution. Adequate possibilities def- 
initely include: a) continuation of A with the help of an 
equation (e.g. demanding the resulting B-field to be po- 
tential or force- free); b) interpolation, either locally (in 
a neighbourhood), or globally (through the whole grid). 
Trying several methods, we concluded that 3-D cubic 
spline interpolation is particularly adequate to the prob- 
lem since it has remarkable advantages over other methods 
(e.g. it does not introduce oscillations in-between grid- 
sites, which would strongly influence the values of the 
derivatives, and it well reproduces the derivatives of an- 
alytically prescribed primary fields). The process of eval- 
uating different continuation methods we went through, 
as well as the comparison of spline interpolation to other 
continuation-methods are described in App. B. 

The 3-D interpolation is performed as three subse- 
quent 1-D interpolations in the three spatial directions 
(Press et al. 1992). For the 1-D splines, we assume nat- 
ural boundaries (the second derivatives are zero at the 
boundaries). Moreover, since in the CA model of LH91 it 
is assumed that around the grid there is a zero field which 
is held constant (see Sec. 2.1), we enlarge the grid by one 
grid point in all directions to include this constant zero- 
layer explicitly, using it however only for the interpolation. 
In the interpolation, the derivatives at the grid-points are 
immediately given by the analytically differentiated inter- 
polating polynomials. 

With the help of this interpolation, the magnetic field 
B and the current J are calculated as derivatives of A, 
according to the MHD prescription: 

B = VAA, (7) 

J = V A B. (8) 

To determine the electric field E, we make the assump- 
tion that imdcr coronal conditions the MHD approach is in 
general valid, and that E is reasonably well approximated 
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by Ohm's law in its simple form, E = rjj —^v/\B, with -q 
the diffusivity and v the fluid velocity. Since the classical 
CA models use no velocity-field, our set-up can yield only 
the resistive part, 

E = r]J. (9) 

In applications such as to solar flares, where the interest 
is in current dissipation events, i.e. in events where rj and 
J are strongly increased, Eq. (9) can be expected to be 
a good approximation to the electric field. Theoretically, 
the convective term in Ohm's law would in general yield 
just a low-intensity, background electric field. 

Eq. (9) needs to be supplemented with a specification 
of the diffusivity rj: Isliker et al. (1998) have shown that 
in the classical CA models the diffusivity adopts the val- 
ues 7/ = 1 at the unstable (bursting) sites, and r] = 
everywhere else. This specifies Eq. (9) completely. 

Rem,ark 1: It is worthwhile noting that, since spline in- 
terpolation has the property to be the least curved of all 
twice differentiable interpolating functions, the grid-size is 
a typical smallest-possible length-scale of field structures, 
or, if we would think the CA to model MHD turbulence, 
it is something like an average smallest possible eddy-size. 

Remark 2: With our set-up, field lines are made available: 

the interpolation was introduced to define derivatives at 
the grid-sites, but it can as well be used to determine the 
vector-potential, and therewith the magnetic field and the 
current, in between the grid sites. However, it is important 
to note that this field in between the grid sites is not used 
for the time-evolution of the model, it merely allows to 
visualize the evolution in the standard way through field- 
lines, if wished so (second order derivatives in-between the 
grid-sites, if needed, would have to be done numerically, 
since else their numerical values would depend on the or- 
der in which the three f D interpolations are done). 

Remark 3: In this paper, we do not change the rules of 
the classical CA models to which we apply our set-up — 
except for the definition of energy release (Sec. 3.1.3). Our 
aim here is to show the usefulness of the set-up and to give 
some results and reveal some aspects by extending pub- 
lished, classical CA models. While the redistribution rules 
have in detail been shown to represent diffusion events 
and fit nicely into an MHD scenario (Isliker et al. 1998), 
the loading process is strongly simplified and poorly fol- 
lows a reasonable flare scenario: For instance, the loading 
process acts independently everywhere in the simulation 
box, whereas according to a realistic flare scenario (see e.g. 
Parker 1993) disturbances should appear independently 
only on one boundary (the photosphere, due to random 
foot-point motions or newly appearing flux), and propa- 
gate then into the interior of the simulation box along the 
magnetic field lines. 

To translate such a realistic loading scenario into the 
language of CA models has not been undertaken, so-far. 



We just note that it would be quite straightforward to in- 
troduce a velocity field into the CA models: e.g. Isliker et 
al. (2000a) propose a CA model which uses a velocity field 
for the loading phase, but this model does not fall into the 
category of classical C A models since it does not follow the 
sand-pile analogy and uses different, MHD based, evolu- 
tion rules. We leave the problem of introducing a velocity 
field and a more physical loading process into the classi- 
cal CA models for a future study. In Isliker et al. 2000b, 
we will — among others — analyze in details what this 
simplified loading process physically represents. 

3. Applications of our set-up 

3.1. Application to the CA model of Lu & Hamilton 
(1991) 

Our first application is to the CA model of LH91 (see 
Sec. 2.1). The LH91 model has a fairly long transient 
phase and reaches finally a stationary state, the so-called 
SOC (self-organized criticality) state, in which spatially 
spreading series of bursts (avalanches) appear, alternating 
with quiet loading phases. The LH91 model gives basically 
three results concerning fiare statistics, namely the distri- 
butions of total energy, peak-flux and durations, which are 
all power-laws with slopes that are in good agreement with 
the observations (Lu et al. 1993, Bromund et al. 1995). 

Superimposing our MHD-frame onto the LH91 model 
such as it stands does not change anyone of the three re- 
sults, since at this first stage we are not interfering with 
the dynamics (i.e. the evolution rules). The set-up allows, 
however, to address several questions in MHD language: 
Our main aim in the subsequent applications is to demon- 
strate that the set-up indeed yields a new and consistent 
interpretation of CA-models, to illustrate the behaviour 
of the secondary variables (currents, magnetic fields), and 
to reveal major features of them. (In the subsequent runs, 
we use a grid of size 30 x 30 x 30, as LH91 did to derive 
their main results.) 

3.1.1. Global structures of the vector-fields 

First, we turn to the question what the global fields 
(vector-potential, magnetic field, current) look like dur- 
ing the SOC state. Thereto, the temporal evolution of 
the model is stopped at an arbitrary time during SOC 
state (in a phase where there are no bursts, i.e. during 
loading), and the magnitude of the fields at a cut with 
fixed z-coordinatc are shown as a function of the x- and 
y-coordinates in Fig. 1. | A(a;, y,z = zo)\ obviously exhibits 
a large-scale organization over the whole grid, it forms a 
global convex surface (Fig. 1(a)). This convex surface has 
a slight random distortion over-lying, which visually can- 
not be discerned in Fig. 1(a), but becomes visible in the 
plot of |B| (Fig. 1(b)), the curl of A, which still exhibits 
large-scale organization all over the grid, but is clearly 
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Fig. 1. Surface and contour plots of the magnitudes of (a) 
the vector-potential {\A\), (b) the magnetic field (|-B|), 
and (c) the current (|J|) as a function of x and y, for 
z = 15 fixed. 

wiggled. Finally, \J\ shows no left-overs of a large scale 
organization anymore, it reflects the random disturbances 
of the convexity of |A| (Fig. 1(c)). 

The large-scale structures shown in Fig. 1 are always 
maintained during the SOC state, neither loading nor 



bursting (and avalanches) destroy them, they just 'trem- 
ble' a little when such events occur. SOC state in the 
extended LH91 model thus implies large-scale organiza- 
tion of the vector-potential and the magnetic field, in the 
characteristic form of Fig. 1. 

The large-scale organization of A is not an artificial 
result of our superimposed set-up, but already inherent in 
the classical LH91 model: in the classical LH91 CA model, 
there is only one variable, the one we call here A, whose 
values arc not affected by the interpolation wc perform 
since it is the primary grid variable, so that Fig. 1(a) is 
true also for the classical, non-extended LH91 model. 

The large scale structure for the primary grid- variable 
I A I is the result of a combined effect: The preferred direc- 
tionality of the loading increments (see Sec. 2.1) tries to in- 
crease \A\ throughout the grid. The redistribution events, 
which already in Bak et al. (1987; 1988) were termed dif- 
fusive events, and which in Isliker et al. (1998) were ana- 
lytically shown to represent local, one-timc-stcp diffusion 
processes, smooth out any too strong spatial unevenness 
of A, and they root the A-field down to the zero level at 
the open boundaries. The result is the convex surface of 
Fig. 1(a), blown- up from below through loading, tied to 
the zero-level at the edges, and forced to a maximum cur- 
vature which is limited by the local, threshold dependent 
diffusion events. 

As the SOC state, so is the large-scale structure of |A| 
independent of the concrete kind of loading, provided it 
fulfills the conditions that the loading increments exhibit 
a preferred directionality and are much smaller than the 
threshold (with symmetric loading, the SOC state is ac- 
tually never reached, see LH91 and Lu et al. (1993)). 

To make sure of the importance of the boundaries, we 
performed runs of the model with closed boundaries, and 
we found that neither a large-scale structure was devel- 
oped in \A\, nor the SOC state was reached. 

3.1.2. Bursts 

To illustrate the role of the current at unstable sites and 
during bursts, we plot in Fig. 2 the magnitude of the cur- 
rent before and after a typical burst: obviously, the current 
at the burst site (x, y, z) = (20, 18, 3) has high intensity 
before the burst (Fig. 2(a)), and is relaxed after the burst 
(Fig. 2(b)). Inspecting a number of other bursts, we found 
that, generally, at sites where the LH91 instability crite- 
rion is fulfilled, the current is increased, too, and that 
bursts dissipate the currents. This is a first hint that clas- 
sical CA models can be interpreted as models for energy 
release through current-dissipation. 

After the burst, at the neighbouring site (21,18,3), 
the intensity of the current is increased, and indeed the 
presented burst gives rise to subsequent bursts, it is one 
event during an avalanche. 

The magnetic field at the bursting site is reshaped, 
in a way which is difficult to interpret when using only 
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Fig. 2. The magnitude of the current (|J|, contour-plot, 
with the ticks pointing 'downhill') as a function of x and 
y at a zoomed cut 2 = 3 through the grid, before (a) and 
after (b) a burst, which occurs in the middle of the plot, 
at {x,y,z) = (20,18,3). 



the magnitude of it {\B\) for visualization. May-be field 
line plots would help visualization, but we leave this for a 
future study. 

3.1.3. Energy release and Ohmic dissipation 

We now turn to the question what relation the energy re- 
lease formula of LH91 (Eq. 6) has to the respective MHD 
relations: In parallel to using the formula of LH91, we de- 
termine the released energy in the following ways, closer to 
MHD: First, we assume it to be proportional to rjJ^ (with 
the diffusivity 77 = 1 at unstable sites, see Sec. 2.2), which 
we linearly interpolate between the two states before and 
after the burst. This is done in two ways, (i) summing over 
the local neighbourhood. 



burst 



t+1 i+1 



dt 
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' ' 1 1 


' ' ' ' 1 










a, 

r; 10-2^ 










< < < < 1 


< < 1 1 
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%eak [arliilrary units] 

Fig. 3. The distribution of the total energy (a), and of the 
peak- flux (b), for different ways of measuring the released 
energy in a burst: Eburst = / dtj^n.n. -^i^iik)^ (solid); 
Eburst — J dt J{xijkY (dotted); Eburst as the difference 
in En n B{xijk,t)'^ before and after the burst (dashed); 
Eburst = jdA^ (dash-dotted). (The distributions are nor- 
malized probability distributions, the last two were shifted 
in both directions for viewing them together with the first 
two.) 



/ ^ <-) n.n.; before ^ *^ n.n.; after ) 



(10) 



and (ii) without summing, but just taking into account 
the current at the central point, 



Hurst = J dtr]J{Xijk,tf 



Ct2 , t2 \ 

V ijk; before ~ ijk; after J 



(11) 



and finally, we monitor the change in magnetic energy due 
to a burst using the difference in magnetic energy in the 
local neighbourhood, 

TP 



hurst 



= E 



{B^^-f°--){Xr,.n.)y 



t n.n. 
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(In Eqs. (10), (11), (12), we assume Ah = 1 and At = 1 for 
the grid-spacing Ah and the time-step At, since, according 
to Isliker et al. (1998), in the classical CA models both 
values are not specified and set to one.) 

The corresponding distributions of total energy and 
peak-flux are shown in Fig. 3, together with the distri- 
butions yielded by the energy-release formula of LII91, 
Eq. (6) (the duration distribution remains the same as in 
the classical LH91 model, namely a power-law, and is not 
shown). Obviously, the four ways of defining the released 
energy give basically similar results, with larger deviations 
only at the low and high energy ends (note that the en- 
ergy in Fig. 3 is in arbitrary units). Using the formula of 
Ohmic dissipation does thus not change the results of the 
classical LH91 model. 

With an estimate of the numerical value of the anoma- 
lous resistivity and of the typical size of a diffusive region 
or the typical difi^usive time, it would be possible to intro- 
duce physical units. We did not undertake this, since all 
three parameters are still known only with large observa- 
tional and theoretical uncertainties. 

3.1.4. The relation of dA to J 

From the similarity of the distributions of the extended 
model with the ones of the classical LH91 model (Fig. 3), 

and from Fig. 2, where it was seen that an instability is 
accompanied by an enhanced current, we are led to ask 
directly for the relation of dA to J, which we plot as a 
function of each other in Fig. 4. Obviously, the two quanti- 
ties are related to each other: above \dA\ « 2, the current 
is an approximate linear function of the stress, around 
\dA\ « 2 the current is zero, and below there is again 
an approximate linear relationship, with negative slope, 
however (above \dA\ « 2 the current J is actually prefer- 
ably along (1,1,1), whereas below it is preferably along 
(— 1, —1, —1), i.e. J is an approximatly linear function of 
dA in the whole range, it merely changes its directivity at 
\dA\ « 2 with respect to dA). In Appendix C, wo show an- 
alytically why with our set-up a more or less close relation 
between dA and J has to be expected. 

Of particular interest in Fig. 4 is that if \dA\ is above 
the threshold Acr = 7, then | J| is also reaching high val- 
ues: obviously, large values of \dA\ imply large values of 
l^j. This confirms the statement made above: The ex- 
tended CA models can be considered as models for en- 
ergy release through current dissipation. It also explains 
why the energy distributions remain very similar when 
the LH91 formula for the amount of energy released in 
a burst (~ dA^, Eq. (6)) is replaced by Ohmic dissipa- 
tion (~ J^, Sec. 3.1.3): bursts occur only for large stresses 
\dA\, where \J\ is also large and an approximate linear 
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Fig. 4. Plot of the magnitude of the current \Jijk\ vs. 
the LII91 stress measure jdAyfel, using the corresponding 
values in the whole grid at a time fixed in a loading phase 
in the SOC state, together with the values at bursting 
sites during an avalanche. 

function of \dA\, so that the distributions of dA^ and 
can be expected to be the same in shape. 

3.2. Application to loading with power-law increments 

Georgoulis and Vlahos (1996, 1998) introduced power-law 
distributed increments for the loading. The main result of 
such a way of driving the system is that the power-law 
indices of the energy-distributions depend on the power- 
law index of the distribution of the loading increments, 
explaining thus the observed variability of the indices 
through the variability of the intensity of the driving. We 
generalize their way of power-law loading, which is for a 
scalar primary field, to a vector field in the following way: 
The anisotropic directivity of the loading increment SA is 
kept (see Sec. 2.1), but \6A\ is now distributed according 
to 

p{\6A\) = C|<5A|-^ (13) 

with \6A\ e [0.01, oo] and /3, the power-law index, a free 

parameter. Simulations were performed for (3 ~ 1.8 and 
(3 = 2.3. Interested in global features implied by the CA 
model, our concern here is the structure of the magnetic 
field. It turns out that the magnetic field exhibits still a 
large scale organization, which is very similar to the one 
of the B-field of the (extended with our set-up) LH91 
model (Fig. 1(b)): for /? = 1.8, the respective plots are 
visually indiscernible, and for /? = 2.3 the overall shape is 
still roughly the same, it merely seems slightly more dis- 
torted. Thus, though the statistical results depend on /?, 
the strength and variability of the loading, the structure 
of the magnetic field remains approximately the same as 
in the case of the extended model of LH91. Large-scale or- 
ganization (in the characteristic form of Fig. 1) must con- 
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for the central point and 




= 4(^) 4- - A 



dA 



(x) 

ijk;n.n. 



Z-^n.n. ijk;n.n. 



(17) 



for those nearest neighbours which fulfill the instabil- 
ity criterion (Eq. 15), where the primed sum is over 
those neighbours for which Eq. (15) holds. The rules for 
j\{y) ^ ^4(2) are completely analogous (so that actually there 

are 18 possibilities to exceed the threshold (Eq. 15) at a 
given site). The released energy is assumed to amount to 



rp{aniso) \ ^ 

rel 



(A^^'> - -A ~)2 



(18) 



Fig. 5. Surface and contour plot of the magnitude of the 
the magnetic field {\B\) as a function of x and y at a 
cut z = 15 through the grid, for the case of anisotropic 
redistribution rules. 



sequently be considered as an inherent property of SOC 
state, through the mechanism explained in Sec. 3.1.1. 

3.3. Application to anisotropic bursts 

Vlahos et al. (1995) introduced anisotropic bursts for solar 
flare CA models, which lead only to small events, but yield 
a steep distribution at small energies, predicting thus a sig- 
nificant over-abundance of small events with a significant 
contribution to coronal heating. We have first to general- 
ize the anisotropic evolution rules, which are again for a 
scalar primary field, to the case of a primary vector field. 
A natural generalization would be to apply the anisotropic 
rules to the absolute magnitude of A, but it turns out that 
this causes the algorithm to get trapped in infinite loops 
(two neighbouring grid-sites trigger each other mutually 
for ever). The same holds if we apply the anisotropic rules 
to the absolute magnitudes of the three components of 
A independently. We finally applied the anisotropic rules 
to the three components of A directly, not using abso- 
lute magnitudes, as also Vlahos et al. (1995) did not use 
absolute magnitudes, and this turned out to lead to a sta- 
tionary asymptotic state: The anisotropic stress in the x- 
component is thus defined as 



dA) 



ijk;n.n. 



:= Alfl - Atl 



ijk 



(14) 



where n.n. stands for one of the six nearest neighbours. 
The instability criterion is 



dA) 



(x) 



> Ac 



ijk;n.n. 

and the redistribution rules become 

.(x) _ .(x) _ 6 



(15) 



(16) 



We performed a run where only the anisotropic burst- 
rules were applied, in order to isolate their effect, although 
the anisotropic burst-rules are used always together with 
the isotropic ones by Vlahos et al. (1995), since alone they 
cannot explain the complete energy distributions of flares. 
In Fig. 5, the magnitude of the magnetic field at a cut 
through the grid is shown (fixed z), for an arbitrary time 
(in the loading phase) during the asymptotic stationary 
state of the model. Clearly, there is no overall large scale 
structure anymore, except that the magnetic field along 
the boundaries is increased. The magnetic field topology 
is thus nearer to the concept of a random, relatively un- 
structured magnetic field than the magnetic field topology 
yielded by the isotropic models in SOC state. 

The anisotropic burst rules do not yield large-scale 
structures, as they are, when used alone, also not able 
to lead the system to SOC state: this is obvious from the 
energy distributions they yield, which are much smaller in 
extent than the ones given by the isotropic rules (see Vla- 
hos et al. 1995), and confirmed by the result of Lu et al. 
(1993) that isotropy of the redistribution rules — at least 
on the average — is a prerequisite to reach SOC state, at 
all. The anisotropic bursts occur independently in all di- 
rections and are in this way not able to organize the field 
in a neighbourhood systematically, and, as a consequence, 
also not in the entire grid. 

The inquiry of the relation of the energy release for- 
mula Eq. (18), which is different from the isotropic formula 
(Eq. 6), to MHD based formulae we leave for a future 
study. We just note that the distributions the anisotropic 
model in our vector-field version yields are at lower ener- 
gies, smaller in extent, and steeper than the ones of the 
isotropic models. 

4. Summary and Conclusions 

Summary 

We have introduced a new set-up for classical solar flare 
CA models which yields, among others, consistency with 

Maxwell's equations (e.g. divergence- free magnetic field), 
and availability of secondary variables such as currents 
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and electric fields in accordance with MHD. Both are new 
for solar flare CA models. The set-up specifies the so far 
open physical interpretation of the CA models. This spec- 
ification is to some extent unavoidably arbitrary, and it 
would definitely be interesting to see what alternative in- 
terpretations would yield — if they can be derived con- 
sistently. We can claim, however, that the interpretation 
we chose is reasonable, it is well-behaved in the sense that 
the derivatives of analytically prescribed vector-potentials 
are reproduced and that the abstract stress-measure of the 
CA models is related to the current, due to general prop- 
erties of spline interpolation. The central problem which 
was to solve is how to calculate derivatives in a CA model, 
i.e. how to continue the primary grid-variable in-between 
the grid sites, since the notion of derivatives is alien in the 
context of CA models quite in general. 

In this article, our main aim with the introduced set- 
up was to demonstrate that the set-up truly extends the 
classical CA models and makes them richer in the sense 
that they contain much more information, now. The main 
features we revealed about the CA models, extended with 
our set-up, are: 

1. Lairge-scale organization of the vector-potential 
and the magnetic field: The field topology during 
SOC state is bound to characteristic large-scale structures 

which span the whole grid, very pronounced for the pri- 
mary grid variable, the vector-potential, but also for the 
magnetic field. Bursts and flares are just slight distur- 
bances propagating over the large-scale structures, which 
are always maintained, also in the largest events. The mag- 
nitude of the current, as a second order derivative of the 
primary field, does not show any obvious large-scale struc- 
ture anymore, it reflects more or less only the random fluc- 
tuations of the large-scale organized magnetic field. It is 
worthwhile noting that the large-scale structure of the pri- 
mary grid- variable is not an artificial result of our set-up, 
but a natural consequence of the SOC state in which the 
system finds itself. The appearance of large-scale struc- 
tures for the primary grid variable was shown here for the 
first time. It may have been known to different authors, 
but it never has explicitly been shown: SOC models for 
fiares are derived in analogy to sand-pile dynamics, and 
the paradigm of a pile reappears in the field topologies of 
the solar flare CA models. 

2. Increased current at unstable grid-sites: Unstable 
sites are characterized by an enhanced current, which is 

reduced after a burst has taken place, as a result of which 
the current at a grid-site in the neighbourhood may be 
increased. 

3. Availability of the electric field: The electric field 

is approximated with the resistive part of Ohm's law in 
its simple form, which can in general be expected to be a 
good approximation in coronal applications and where the 
interest is in current- dissipation events, e.g. in the case of 
solar flares. 



4. Energy release in terms of Ohmic dissipation: 

We replaced the some-what ad hoc formula in the CA 
models to estimate the energy released in a burst with the 
expression for Ohmic dissipation in terms of the current. 
The distributions yielded in this way are very similar to 
the ones based on the ad hoc formula, so that the results 
of the CA models remain basically unchanged. 

5. CA as models for current dissipations: As a con- 
sequence of point 2 and 4 in this list, and of the fact that 

there is an approximate linear relation between the cur- 
rent and the stress measure of the CA, we can conclude 
that the extended CA models can be considered as models 
for energy release through current dissipation. 

4-2. Conclusions 

Our set-up is to be contrasted to the recently suggested 
MHD-derived (not based on the sand-pile analogy) CA 
models of Einaudi & Velli (1999), MacPherson & MacK- 
innon (1999), Longcope and Noonan (2000), and Isliker 
et al. (2000a). They all suggest new evolution rules, de- 
rived from MHD, and all in different ways (they actu- 
ally focus on different processes, namely the microscopic, 
macroscopic, and mesoscopic physics, respectively, in ac- 
tive regions). Our set-up, on the other hand, uses existing 
CA models, does not interfere (if not wished) with their 
evolution rules, does also not change their main results, 
as shown, but reinterprets them, extends them essentially, 
and makes them compatible with MHD. 

The set-up wc introduced allows different future appli- 
cations and posing questions which could not be asked so 
far in the frame of CA models. In preparation is a study 
(Isliker et al. 2000b) to reveal in detail what physical flare 
scenario the extended CA models imply. We will address 
the questions: (1) how to interpret the small scale pro- 
cesses of the models (loading and bursting) in terms of 
MHD; (2) what the global flare-scenario implied by the 
models is; (3) whether the global magnetic field topol- 
ogy of the models can be considered to represent observed 
magnetic topologies in active regions; (4) what spatio- 
temporal evolution of the electric field during flares is 
yielded by the models. 

A different future application we plan with CA models 
extended with our set-up is the introduction of particles 
into the models, with the aim to study thermal emission, 
particle acceleration, and non-thermal emission. This will 
allow a much deeper comparison of the CA models to ob- 
servations than was possible so far, and this is actually the 
most important beneflt of the set-up we introduced. Such 
comparisons will allow a new judgment of the adequate- 
ness or not of classical CA models (in their current form) 
to the problem of solar flares, beyond the three statistical 
distributions of the primarily released energy. Solar flare 
CA models which include particle acceleration would rep- 
resent the first global and complete model for solar flares. 
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Appendix A: Temporal evolution of the CA 

The temporal evolution of the CA models presented in 
this article is governed by the following rules: 

0. initializing 

1. loading 

2. scanning: create a list of the unstable sites; if there are 
none, return to loading (1) 

3. scanning and bursting: redistribute the fields at the 
unstable sites which are in the list created in the scan- 
nings 2 or 4 

4. scanning: create a list of the unstable sites. If there are 
any, go to bursting (3), else return to loading (1) 

The extra scannings 2 and 4 are needed for causality: if a 
site becomes unstable through a burst in the neighbour- 
hood, then it should be redistributed in the subsequent 
scan, and not in the same as the primary unstable site. 
The same is true for the scanning 4, since in the next 
bursting phase (if any) only those sites should burst who 
had become unstable through a burst in their neighbour- 
hood during the foregoing time-step. 

As a time-step is considered one scanning of the grid, 
point 3. The released energy per time-step is the sum of 
all the energy released by bursts in this time-step (a burst 
is considered a single redistribution event in 3). We term 
a flare or avalanche the loop 3,4, from the occurring of 
the first burst in 3 until the activity has died out and one 
returns via the scanning 4 to loading (1). The duration of 
the flare is the number of time-steps it lasted, the total 
flare energy is the sum of all the energies released in the 
duration of the flare, and the peak flux or peak energy is 
the maximum of the energies of all the time-steps of the 
flare. 

Appendix B: Why spline interpolation is 

particularly adequate: comparison 
to other methods of continuation 

We mentioned in Sec. 2.2 that other possibilities for con- 
tinuation of the vector-potential besides spline interpola- 
tion would be: a) continuation of A with the help of an 
equation; b) other kinds of interpolation, either locally (in 
a neighbourhood), or globally (through the whole grid). 
Possibility a) implies that an equation has to be solved in 
each time-step (after each loading and after each burst), in 
the worst case numerically, with open boundary condition 
and the Aijk given at the grid-sites. This computational 



effort might slow down the algorithm of the model con- 
siderably (and bring it near to the computational effort 
of MHD equation integration). Besides that, the problem 
is what equation to use: to make the magnetic field al- 
ways a potential field (i.e. using a corresponding equation 
for the vector-potential A) implies that, from the point 
of view of MHD, at all times a very 'well-behaved' mag- 
netic field resides in the CA, with no tendency towards 
instabilities, which makes it diflacult to understand why 
bursts should occur at all, since critical quantities such as 
currents do not become excited. A better candidate could 
be expected to be force-freeness, except that, possibly, one 
may be confronted with incompatibility of the boundary 
conditions with the vector-potential values given at the 
grid-sites, i.e. existence-problems for solutions eventually 
arise. 

Though definitely possibility a) cannot be ruled out 
on solid grounds, we found it more promising to proceed 
with possibility b), interpolation. A guide-line for choos- 
ing a particular interpolation method is the reasonable 
demand that the interpolation should not introduce wild 
oscillations in-between grid-sites, for wc want to assure 
that the derivatives at the grid sites, which are very sen- 
sitive to such oscillations, are not 'random' values solely 
due to the interpolation, but that they reflect more or 
less directly the change of the primary grid-variable from 
grid-site to grid-site. This calls for interpolating functions 
which are as little curved as possible. 

The ciasiest and fastest way of interpolating would be 
to perform local interpolations around a point and its 
nearest neighbours (e.g. using low-order polynomials or 
trigonometric functions of different degrees). This inter- 
polation leads, however, to ambiguities for the derivatives: 
the derivatives, say at a point Xijk, are not the same, if 
the used interpolation is centered at Xijk, with the ones 
calculated with an interpolation centered at e.g. Xi^ijk- 
In this sense, local interpolation is not self-consistent, the 
derivatives at a grid-site depend on where the used inter- 
polation is centered. 

Finally, we are left with global interpolation through 
the whole grid. Among the candidates are, besides more 
exotic interpolating functions, polynomials of degree equal 
to the grid size, trigonometric functions (also in the 
form of Fourier-transforms) , low-order smooth polynomi- 
als (e.g. splines). The first candidate, polynomials of a 
high degree n (with n the number of grid points in one di- 
rection), we reject immediately since it is notorious for its 
strong oscillations in-between grid-sites, mainly towards 
the edges of the grid. We tried the second candidate, 
trigonometric interpolation, in the form of discrete Fourier 
transform. Testing this by prescribing analytic functions 
for A(x) and comparing the numerical derivatives with 
the analytic ones, it turned out that there arise problems 
with representing structures in A as large as the entire 
grid (the wave- number spectrum is too limited), and with 
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structures as short as roughly the grid-spacing (different 
prescribed short structures are taken for the same). 

Trying cubic spline- interpolation, we found that it does 
not suffer from the problems stated for the other types of 
interpolation: neither does it introduce wild oscillations, 
unmotivated by the values at the grid-sites, nor does spline 
interpolation have problems with describing large or small 
scale structures (if a functional form of A is prescribed, 
then the analytic derivatives and the derivatives yielded 
by the interpolation give very close values, in general). 

Moreover, based on results of Sec. 3, App. C, and 
Isliker et al. (1998), there is another reason why spline- 
interpolation is particularly adequate to our problem: It 
relates the quantity dA (Eq. (2)), which measures the 
stress at a site in the CA model, closely toS/^A, the Lapla- 
cian of A (see App. C). The latter is related to the cur- 
rent (J = -^V^A ^V(VA)), which, from the point 
of view of MHD, can be considered as a measure of stress 
in the magnetic field configuration. If this relation would 
not hold, then the redistribution rules (Eqs. (4) and (5)) 
of the CA would not be interpretable as the diffusion pro- 
cess revealed by Isliker et al. (1998), and the instability 
criterion (Eq. 3) would not be so closely related to the 
current (see Sec. 3 and App. C). 

B.l. Why in particular differencing is not adequate to 

calculate derivatives in a CA 

We had rejected above (Sec. 1, Sec. 2.2) the use of dif- 
ference expressions to calculate derivatives, stating that 
differencing is not in the spirit of CA models quite in gen- 
eral, since the nature of CA is truly discrete. We think it 
worthwhile to make this argument more concrete and to 
show what problems arise if differencing were used: 

1. Consistency with the evolution rules: Isliker et al. (1998) 
have shown that the classical solar flare CA are not just 
the discretized form of a differential equations. Instead, 
they describe the time-evolution of a system by rules 
which express the direct transition from a given initial to 
a final state which is the asymptotic solution of a simple 
diffusion equation. The time-step corresponds therewith 
to the average time needed for smallest scale structures 
(structures as large as a neighbourhood) to diffuse, and 
the grid-size corresponds to the size of these smallest oc- 
curring structures. Assuming that the CA models were 
just discretized differential equations would lead to severe 
mathematical and physical contradictions and inconsis- 
tencies (continuity for Ah ^ is violated (with Ah the 
grid-size), and negative diffusivities appear). Therewith, 
in order to be consistent with the evolution rules, which 
assume a finite grid-size, one cannot assume for the pur- 
pose of differentiating this same grid-size to be approxi- 
mately infinitesimal. 

2. Derivatives as difference expressions are not self- 
consistent: There are several equivalent ways to define 
numerical derivatives with the use of difference expres- 



sions: there are e.g. the backward difference dxAx(xijk) = 
{Ax{xijk) — Ax{xi-ijk)) / Ah, and the forward difference 
da:Ax{xijk) = {Ax{xi+ijk) - Ax{xijk))/Ah. Both should 
give comparable values in a given application, else, in 
the context of differential equation integration, one would 
have to make the resolution higher. In the case of CA- 
modcls, we find that the two difference expressions yield 
values which differ substantially from each other: E.g. for 
an initial loading of the grid with independent random val- 
ues for the j4-field, the difference between the backward 
and the forward difference expression can be as large as 
the field itself. Such an initial condition would of course 
not make sense in the context of partial differential equa- 
tions, in the context of CA, however, it is a reasonable 
starting configuration, and the evolution is unaffected by 
such an initial loading. Moreover, when the CA models we 
discuss in this article have reached the SOC state, then 
the differences between e.g. the backward- and forward- 
difference expressions can be as large as 400%. There is 
no way to reduce this discrepancy, since grid-refinement 
is principally impossible for CA: the evolution is governed 
by a set of rules, and making the grid spacing smaller by 
introducing new grid-points in-between the old ones would 
actually just mean to make the grid larger, since the evo- 
lution rules remain the same, there are no rules for half 
the grid-spacing. 

Appendix C: Relation of dA to AA 

The stress measure of LH91, dAijk = Aijk — 
^ ri -^n.n., Can bc related to continuous expressions 
by representing the values of as Taylor-series ex- 

pansions around Xijk, setting the spatial differences to 
Ah = 1. It turns out that e.g. 

dA^ = -IaA, - ^{dt + d^y + dt)A, - ... (C.l) 

and so on for the other two components. In general, it 
is therefore not adequate to consider dA to be a good 
4th order approximation to A^l, since higher order cor- 
rections can be large, they depend on the way the vector 
potential is continued in-between grid-sites. If we had, for 
instance, chosen global polynomial interpolation instead 
of spline-interpolation, the higher order terms would not 
be negligible, above all towards the edges of the grid, since 
polynomial interpolation is known for introducing fluctu- 
ations near the edges of the grid. Consequently, dA would 
be a bad approximation to A A. In order dA to be a good 
approximation to AA, interpolation with, for example, 
3rd order polynomials would be an optimum choice {dA 
would be an exact approximation to AA). Thus, 3rd order 
polynomials would be the choice for local interpolation, 
which, however, is not applicable, since it introduces dis- 
continuities in B and J (see App. B). The way out of the 
dilemma we suggested in this article is the use of cubic 
splines, which provide global interpolation with 3rd order 
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polynomials, with B and J continuous, and only third or- 
der derivatives are discontinuous (this is the price of the 
compromise). For splines then, Eq. (C.l) writes as 



Vlahos, L., Georgoulis, M., Kluiving, R., Paschos, P., 1995, A&A 
299, 897 



( d^At-dlA;) 
( d^At-d^A-) 



(C.2) 



due to the discontinuities in the 3rd order derivatives (the 
superscripts + and — refer to the right and left derivative, 
respectively). Thus, in case where the third order right 
and left derivatives are not too different, dA is a good 
approximation to AA. 
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